
ValueAdd <- pcp_claims[, .(FirstYear=min(as.numeric(VisitYear)), LastYear=max(as.numeric(VisitYear)), FrequentYear=as.numeric(names(sort(table(VisitYear), dec=T)[1])), N=length(VisitSID),
                           ProviderType=names(sort(table(ProviderType), dec=T)[1]), Sta3n=names(sort(table(Sta3n), dec=T)[1]), InstitutionSID=names(sort(table(InstitutionSID), dec=T)[1])), by="StaffSSN"]

ValueAdd$StaffSSN <- as.numeric(as.character(ValueAdd$StaffSSN))

FE <- drift[, .(FE_MH_bare=mean(FE_MH_bare), FE_MH=mean(FE_MH), FE_PH_bare=mean(FE_PH_bare), FE_PH=mean(FE_PH), FE_ACSC_bare=mean(FE_ACSC_bare), FE_ACSC=mean(FE_ACSC), AnnualCaseCount=length(PatientICN)), by="StaffSSN"]
FE$StaffSSN <- as.numeric(as.character(FE$StaffSSN))
ValueAdd <- merge(x=ValueAdd, y=FE, by.x="StaffSSN", by.y="StaffSSN")

fwrite(ValueAdd, "~/ValueAdd.csv")

###############################################################
#######         FIGURE 1 HISTOGRAM OF VALUE ADD         #######
###############################################################

par(mfrow=c(3,1), mai=c(0.5,0.9,0.7,0.2))
ValueAdd$FE_MH <- ValueAdd$FE_MH-mean(ValueAdd$FE_MH)
ValueAdd$FE_PH <- ValueAdd$FE_PH-mean(ValueAdd$FE_PH)
ValueAdd$FE_ACSC <- ValueAdd$FE_ACSC-mean(ValueAdd$FE_ACSC)
hist(-ValueAdd$FE_MH[which(abs(ValueAdd$FE_MH)<=0.062)], cex.axis=2, cex.lab=2, cex.main=2, main="Mental Health", xlab="", xlim=c(-0.062,0.062), breaks=seq(-0.062,0.062,length.out=40), col=alpha("royalblue4",0.5))
hist(-ValueAdd$FE_PH[which(abs(ValueAdd$FE_PH)<=0.062)], cex.axis=2, cex.lab=2, cex.main=2, main="Circulatory Conditions", xlab="", xlim=c(-0.062,0.062), breaks=seq(-0.062,0.062,length.out=40), col=alpha("firebrick3",0.5))
hist(-ValueAdd$FE_ACSC[which(abs(ValueAdd$FE_ACSC)<=0.062)], cex.axis=2, cex.lab=2, cex.main=2, main="ACSC", xlab="", breaks=seq(-0.062,0.062,length.out=40), xlim=c(-0.062,0.062), col=alpha("darkorange",0.5))

# Standardize and take negative
ValueAdd$FE_MH_std <- -scale(ValueAdd$FE_MH)
ValueAdd$FE_MH_bare_std <- -scale(ValueAdd$FE_MH_bare)
ValueAdd$FE_PH_std <- -scale(ValueAdd$FE_PH)
ValueAdd$FE_PH_bare_std <- -scale(ValueAdd$FE_PH_bare)
ValueAdd$FE_ACSC_std <- -scale(ValueAdd$FE_ACSC)
ValueAdd$FE_ACSC_bare_std <- -scale(ValueAdd$FE_ACSC_bare)



##############################################################################
############           TABLE 1: SUMMARY STATISTICS               #############
##############################################################################

pcp_claims$relationshipLength[which(is.na(pcp_claims$relationshipLength))] <- as.numeric(as.Date("2020-09-01", "%Y-%m-%d")-as.Date(pcp_claims$RelationshipStartDate[which(is.na(pcp_claims$relationshipLength))]))

pcp_claims$length_truncate <- pcp_claims$relationshipLength
pcp_claims$length_truncate[which(pcp_claims$relationshipLength>=(365*3))] <- 365*3
pcp_claims$length_truncate[which(is.na(pcp_claims$relationshipLength))] <- 365*3

pcp_claims$Black <- (pcp_claims$RaceEthnicity=="Black")
pcp_claims$API <- (pcp_claims$RaceEthnicity=="API")
pcp_claims$White <- (pcp_claims$RaceEthnicity=="White")
pcp_claims$Hispanic <- (pcp_claims$RaceEthnicity=="Hispanic")
pcp_claims$Native <- (pcp_claims$RaceEthnicity=="Native")
pcp_claims$MarriedNow <- (pcp_claims$Married=="1")
pcp_claims$MarriedPrevious <- (pcp_claims$Married=="2")
pcp_claims$MarriedNever <- (pcp_claims$Married=="3")
pcp_claims$Korean <- (pcp_claims$PeriodOfService=="Korean")
pcp_claims$Vietnam <- (pcp_claims$PeriodOfService=="Vietnam")
pcp_claims$Gulf <- (pcp_claims$PeriodOfService=="Gulf")
pcp_claims$OtherPeriod <- (!pcp_claims$PeriodOfService %in% c("Korean", "Vietnam", "Gulf"))
pcp_claims$Muscu <- (pcp_claims$MDC=="Diseases of the musculoskeletal system and connective tissue")
pcp_claims$Circulatory <- (pcp_claims$MDC=="Diseases of the circulatory system")
pcp_claims$Endocrine <- (pcp_claims$MDC=="Endocrine; nutritional; and metabolic diseases and immunity disorders")
pcp_claims$Mental <- (pcp_claims$MDC=="Mental Illness")
pcp_claims$Respiratory <- (pcp_claims$MDC=="Diseases of the respiratory system")
pcp_claims$OtherMDC <- (!pcp_claims$MDC %in% c("Diseases of the musculoskeletal system and connective tissue", "Diseases of the circulatory system", "Endocrine; nutritional; and metabolic diseases and immunity disorders", "Mental Illness", "Diseases of the respiratory system"))
pcp_claims$AnnualVACheckAmmount <- pcp_claims$AnnualVACheckAmmount*100/pcp_claims$cpi

pcp_claims$FE_PH_std_res <- resid(felm(FE_PH_std~0|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins, data=pcp_claims))
bottom_tercile <- quantile(pcp_claims$FE_PH_std_res, 0.33333)
top_tercile <- quantile(pcp_claims$FE_PH_std_res, 0.66666)

bottom_tercile <- ValueAdd$StaffSSN[which(ValueAdd$FE_PH_std<=quantile(ValueAdd$FE_PH_std,0.33))]
middle_tercile <- ValueAdd$StaffSSN[which(ValueAdd$FE_PH_std>quantile(ValueAdd$FE_PH,0.33) & ValueAdd$FE_PH_std<=quantile(ValueAdd$FE_PH_std,0.66))]
top_tercile <- ValueAdd$StaffSSN[which(ValueAdd$FE_PH_std>quantile(ValueAdd$FE_PH_std,0.66))]

round(
  cbind(
    colMeans(pcp_claims[, c("Age", "API", "Black", "Hispanic", "Native", "White", "MarriedNow", "MarriedPrevious", "MarriedNever", "PatientIncome", "Medicare", "Medicaid", "Korean", "Vietnam", "Gulf", "OtherPeriod", "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount",
                            "PriorVisit", "MH_dummy_lag1_v3", "PH_dummy_lag1_v3", "ACSC_lag1", "WaitTime", "Muscu", "Circulatory", "Endocrine", "Mental", "Respiratory", "OtherMDC", "length_truncate")], na.rm=T)*100,
    colMeans(pcp_claims[StaffSSN %in% bottom_tercile, c("Age", "API", "Black", "Hispanic", "Native", "White", "MarriedNow", "MarriedPrevious", "MarriedNever", "PatientIncome", "Medicare", "Medicaid", "Korean", "Vietnam", "Gulf", "OtherPeriod", "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount",
                                                        "PriorVisit", "MH_dummy_lag1_v3", "PH_dummy_lag1_v3", "ACSC_lag1", "WaitTime", "Muscu", "Circulatory", "Endocrine", "Mental", "Respiratory", "OtherMDC", "length_truncate")], na.rm=T)*100,
    colMeans(pcp_claims[StaffSSN %in% middle_tercile, c("Age", "API", "Black", "Hispanic", "Native", "White", "MarriedNow", "MarriedPrevious", "MarriedNever", "PatientIncome","Medicare", "Medicaid", "Korean", "Vietnam", "Gulf", "OtherPeriod", "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount",
                                                        "PriorVisit", "MH_dummy_lag1_v3", "PH_dummy_lag1_v3", "ACSC_lag1", "WaitTime", "Muscu", "Circulatory", "Endocrine", "Mental", "Respiratory", "OtherMDC", "length_truncate")], na.rm=T)*100,
    colMeans(pcp_claims[StaffSSN %in% top_tercile, c("Age", "API", "Black", "Hispanic", "Native", "White", "MarriedNow", "MarriedPrevious", "MarriedNever", "PatientIncome", "Medicare", "Medicaid", "Korean", "Vietnam", "Gulf", "OtherPeriod", "ServiceConnectedFlag", "UnemployableFlag", "ExposedToAgentOrangeFlag", "RadiationExposureIndicatedFlag", "AnnualVACheckAmount",
                                                     "PriorVisit", "MH_dummy_lag1_v3", "PH_dummy_lag1_v3", "ACSC_lag1", "WaitTime", "Muscu", "Circulatory", "Endocrine", "Mental", "Respiratory", "OtherMDC", "length_truncate")], na.rm=T)*100), 1)


###############################################################
#######         FIGURE 2 BALANCE TABLES/FIGURES         #######
###############################################################


elixhauser <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 
DROP TABLE IF EXISTS #ICD9
DROP TABLE IF EXISTS #ICD10
DROP TABLE IF EXISTS #relevant_dx
DROP TABLE IF EXISTS #cohort;

SELECT PatientICN, PatientSID, VisitDateTime, VisitSID, StaffSSN INTO #cohort FROM OMHO_QFR.ECON.appointments


--- DIAGNOSIS CODES
-- including CMS

-- ICD 9
SELECT DISTINCT ICDCode, ICD9SID, type
INTO #ICD9
FROM(
	SELECT 
		a.ICD9Code AS ICDCode, a.ICD9SID,
		CASE
		WHEN a.ICD9Code IN ('E850.1', 'E850.2', 'E935.1', 'E935.2', 'E980.0', '965.00', '965.02', '965.09') THEN 'opioid overdose'

		WHEN a.ICD9Code IN ('E850.0', 'E935.0', '965.01') THEN 'heroin overdose'

		WHEN OUD=1 THEN 'opioid dependence'
	
		WHEN SuicideAttempt=1 THEN 'suicide'

		WHEN DEPRESS=1 THEN 'depression'

		WHEN homeless=1 THEN 'homeless'
		
		WHEN SAE_Falls=1 THEN 'fall'

		WHEN SUD_Active_Dx=1 THEN 'SUD'

		WHEN SAE_Acet=1 OR SAE_Falls=1 OR SAE_OtherAccident=1 OR SAE_OtherDrug=1 OR SAE_sed=1 OR SAE_Vehicle=1 THEN 'SAE'

		-- Elixhauser-related ones
		WHEN a.ICD9Code IN ('398.91', '402.01', '402.11', '402.91', '404.01', '404.03', '404.11', '404.13', '404.91', '404.93') THEN 'congestive heart failure'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('426.0', '426.7', '426.9', '426.10', '427.0', '427.1', '427.2', '427.3', '427.4', '427.6', '427.7', '427.8', '427.9', '785.0') 
		OR a.ICD9Code IN ('426.10', '426.12', '426.13') THEN 'cardiac arrhythmia'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('093.2', '746.3', '746.4', '746.5', '746.6', 'V42.2', 'V43.3')
		OR a.ICD9Code LIKE '394.%' OR a.ICD9Code LIKE '395.%' OR a.ICD9Code LIKE '396.%' OR a.ICD9Code LIKE '397.%' 
		OR a.ICD9Code LIKE '746.3%' OR a.ICD9Code LIKE '746.4%' OR a.ICD9Code LIKE '746.5%' OR a.ICD9Code LIKE '746.6%' THEN 'valvular disease'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('415.0', '415.1', '417.0', '417.8', '417.9') OR a.ICD9Code LIKE '416.%' THEN 'pulmonary circulation disorders'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('093.0', '437.3', '443.1', '443.2', '443.3', '443.4', '443.5', '443.6', '443.7', '443.8', '443.9', '447.1', '557.1', '557.9')
		OR a.ICD9Code LIKE '440.%' OR a.ICD9Code LIKE '441.%' THEN 'peripheral vascular disorders' 
	
		WHEN a.ICD9Code LIKE '401.%' THEN 'hypertension, uncomplicated' 

		WHEN a.ICD9Code LIKE '402.%' OR a.ICD9Code LIKE '403.%' OR a.ICD9Code LIKE '404.%' OR a.ICD9Code LIKE '405.%' THEN 'hypertension, complicated'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('334.1', '344.0', '344.1', '344.2', '344.3', '344.4', '344.5', '344.6', '344.9') 
		OR a.ICD9Code LIKE '342.%' OR a.ICD9Code LIKE '343.%' THEN 'paralysis'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('331.9', '332.0', '332.1', '333.4', '333.5', '336.2', '348.1', '348.3', '780.3', '784.3')
		OR a.ICD9Code='333.92'
		OR a.ICD9Code LIKE '334.%' OR a.ICD9Code LIKE '335.%' OR a.ICD9Code LIKE '340.%' OR a.ICD9Code LIKE '341.%' OR a.ICD9Code LIKE '345.%' THEN 'other neurological disorders'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('416.8', '416.9', '506.4', '508.1', '508.8')
		OR a.ICD9Code LIKE '490.%' OR a.ICD9Code LIKE '491.%' OR a.ICD9Code LIKE '492.%' OR a.ICD9Code LIKE '493.%' OR a.ICD9Code LIKE '494.%' 
		OR a.ICD9Code LIKE '495.%' OR a.ICD9Code LIKE '496.%' OR a.ICD9Code LIKE '497.%' OR a.ICD9Code LIKE '498.%' OR a.ICD9Code LIKE '499.%' 
		OR a.ICD9Code LIKE '500.%' OR a.ICD9Code LIKE '501.%' OR a.ICD9Code LIKE '502.%' OR a.ICD9Code LIKE '503.%' OR a.ICD9Code LIKE '504.%' OR a.ICD9Code LIKE '505.%'THEN 'chronic pulmonary disease'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('250.0', '250.1', '250.2', '250.3') THEN 'diabetes, uncomplicated'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('250.4', '250.5', '250.6', '250.7', '250.8', '250.9') THEN 'diabetes, complicated'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('240.9', '246.1', '246.8') OR a.ICD9Code LIKE '243.%' OR a.ICD9Code LIKE '244.%' THEN 'hypothyroidism'

		WHEN a.ICD9Code IN ('403.01', '403.11', '403.91', '404.01', '404.03', '404.12', '404.13', '404.92', '404.93', '588.0', 'V42.0', 'V45.1')
		OR a.ICD9Code LIKE '585.%' OR a.ICD9Code LIKE '586.%' OR a.ICD9Code LIKE 'V56.%' THEN 'renal failure'

		WHEN a.ICD9Code IN ('070.22', '070.23', '070.32', '070.33', '070.44', '070.54', '070.6', '070.9', '456.0', '456.1', '456.2',
		'572.2', '572.3', '572.4', '572.5', '572.6', '572.7', '572.8', '573.3', '573.4', '573.8', '573.9', 'V42.7')
		OR a.ICD9Code LIKE '570.%' OR a.ICD9Code Like '571.%' THEN 'liver disease'

		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('531.7', '531.9', '532.7', '532.9', '533.7', '533.9', '534.7', '534.9') THEN 'peptic ulcer disease'

		WHEN a.ICD9Code LIKE '042.%' OR a.ICD9Code LIKE '043.%' OR a.ICD9Code LIKE '044.%' THEN 'aids/hiv'

		WHEN a.ICD9Code IN ('203.0', '238.6') OR a.ICD9Code LIKE '200.%' OR a.ICD9Code LIKE '201.%' OR a.ICD9Code LIKE '202.%' THEN 'lymphoma'

		WHEN a.ICD9Code LIKE '196.%' OR a.ICD9Code LIKE '197.%' OR a.ICD9Code LIKE '198.%' OR a.ICD9Code LIKE '199.%' THEN 'metastatic cancer'

		WHEN a.ICD9Code LIKE '140.%' OR a.ICD9Code LIKE '141.%' OR a.ICD9Code LIKE '142.%' OR a.ICD9Code LIKE '143.%' OR a.ICD9Code LIKE '144.%' OR a.ICD9Code LIKE '145.%' OR a.ICD9Code LIKE '146.%' OR a.ICD9Code LIKE '147.%' OR a.ICD9Code LIKE '148.%' OR a.ICD9Code LIKE '149.%' OR 
		a.ICD9Code LIKE '150.%' OR a.ICD9Code LIKE '151.%' OR a.ICD9Code LIKE '152.%' OR a.ICD9Code LIKE '153.%' OR a.ICD9Code LIKE '154.%' OR a.ICD9Code LIKE '155.%' OR a.ICD9Code LIKE '156.%' OR a.ICD9Code LIKE '157.%' OR a.ICD9Code LIKE '158.%' OR a.ICD9Code LIKE '159.%' OR 
		a.ICD9Code LIKE '160.%' OR a.ICD9Code LIKE '161.%' OR a.ICD9Code LIKE '162.%' OR a.ICD9Code LIKE '163.%' OR a.ICD9Code LIKE '164.%' OR a.ICD9Code LIKE '165.%' OR a.ICD9Code LIKE '166.%' OR a.ICD9Code LIKE '167.%' OR a.ICD9Code LIKE '168.%' OR a.ICD9Code LIKE '169.%' OR 
		a.ICD9Code LIKE '170.%' OR a.ICD9Code LIKE '171.%' OR a.ICD9Code LIKE '172.%' OR a.ICD9Code LIKE '174.%' OR a.ICD9Code LIKE '175.%' OR a.ICD9Code LIKE '176.%' OR a.ICD9Code LIKE '177.%' OR a.ICD9Code LIKE '178.%' OR a.ICD9Code LIKE '179.%' OR 
		a.ICD9Code LIKE '180.%' OR a.ICD9Code LIKE '181.%' OR a.ICD9Code LIKE '182.%' OR a.ICD9Code LIKE '183.%' OR a.ICD9Code LIKE '184.%' OR a.ICD9Code LIKE '185.%' OR a.ICD9Code LIKE '186.%' OR a.ICD9Code LIKE '187.%' OR a.ICD9Code LIKE '188.%' OR a.ICD9Code LIKE '189.%' OR 
		a.ICD9Code LIKE '190.%' OR a.ICD9Code LIKE '191.%' OR a.ICD9Code LIKE '192.%' OR a.ICD9Code LIKE '193.%' OR a.ICD9Code LIKE '194.%' OR a.ICD9Code LIKE '195.%' THEN 'solid tumour without metastasis'
	
		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('701.0', '710.0', '710.1', '710.2', '710.3', '710.4', '710.8', '710.9', '711.2', '719.3', '728.5')OR a.ICD9Code IN ('728.89', '729.30')
		OR a.ICD9Code LIKE '446.%' OR a.ICD9Code LIKE '714.%' OR a.ICD9Code LIKE '720.%' OR a.ICD9Code LIKE '725.%' THEN 'rheumatoid arthritis'
	
		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('287.1', '287.3', '287.4', '287.5') OR a.ICD9Code LIKE '286.%' THEN 'coagulopathy'

		WHEN a.ICD9Code='278.0' THEN 'obesity'
	
		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('783.2', '799.4') OR a.ICD9Code LIKE '260.%' OR a.ICD9Code LIKE '261.%' OR a.ICD9Code LIKE '262.%' OR a.ICD9Code LIKE '263.%' THEN 'weight loss'
	
		WHEN SUBSTRING(a.ICD9Code, 1, 5)='253.6' OR a.ICD9Code LIKE '276.%' THEN 'fluid and electrolyte disorders'

		WHEN SUBSTRING(a.ICD9Code, 1, 5)='280.0' THEN 'blood loss anaemia'
	
		WHEN SUBSTRING(a.ICD9Code, 1, 5) IN ('280.1', '280.2', '280.3', '280.4', '280.5', '280.6', '280.7', '280.8', '280.9') OR a.ICD9Code LIKE '281.%' THEN 'deficiency anaemia'

		ELSE NULL END AS type
	FROM CDWWork.Dim.ICD9 AS a
	LEFT OUTER JOIN OMHSP_PERC_Share.DOEx.Lookup_ICD9 AS lookup
		ON a.ICD9SID=lookup.ICD9SID
) AS z
WHERE type IS NOT NULL


----- ICD 10
SELECT DISTINCT ICDCode, ICD10SID, type
INTO #ICD10
FROM(
	SELECT 
		a.ICD10Code AS ICDCode, a.ICD10SID,
		CASE
		WHEN a.ICD10Code LIKE 'T40.0X%' OR a.ICD10Code LIKE 'T40.2X%' OR a.ICD10Code LIKE 'T40.3X%'  OR a.ICD10Code LIKE 'T40.4X%' OR a.ICD10Code LIKE 'T40.6%' THEN 'opioid overdose'

		WHEN a.ICD10Code LIKE 'T40.1X%' THEN 'heroin overdose'

		WHEN OUD=1 THEN 'opioid dependence'
	
		WHEN SuicideAttempt=1 THEN 'suicide'

		WHEN DEPRESS=1 THEN 'depression'

		WHEN homeless=1 THEN 'homeless'
		
		WHEN SAE_Falls=1 THEN 'fall'

		WHEN SUD_Active_Dx=1 THEN 'SUD'

		WHEN SAE_Acet=1 OR SAE_Falls=1 OR SAE_OtherAccident=1 OR SAE_OtherDrug=1 OR SAE_sed=1 OR SAE_Vehicle=1 THEN 'SAE'

		-- Elixhauser-related ones
		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I09.9', 'I11.0', 'I13.0', 'I13.2', 'I25.5', 'I42.0', 'I42.5', 'I42.6', 'I42.7', 'I42.8', 'I42.9', 'P29.0')
		OR a.ICD10Code LIKE 'I43.%' OR a.ICD10Code LIKE 'I50.%' THEN 'congestive heart failure'
		
		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I44.1', 'I44.2', 'I44.3', 'I45.6', 'I45.9', 'R00.0', 'R00.1', 'R00.8', 'T82.1', 'Z45.0', 'Z95.0')
		OR a.ICD10Code LIKE 'I47.%' OR a.ICD10Code LIKE 'I48.%' OR a.ICD10Code LIKE 'I49.%' THEN 'cardiac arrhythmia' 

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('A52.0', 'I09.1', 'Q23.0', 'Q23.1', 'Q23.2', 'Q23.3', 'Z95.2', 'Z95.3', 'Z95.4')
		OR a.ICD10Code LIKE 'I05.%' OR a.ICD10Code LIKE 'I06.%' OR a.ICD10Code LIKE 'I07.%' OR a.ICD10Code LIKE 'I08.%' OR a.ICD10Code LIKE 'I34.%' OR a.ICD10Code LIKE 'I35.%' OR a.ICD10Code LIKE 'I36.%' 
		OR a.ICD10Code LIKE 'I37.%' OR a.ICD10Code LIKE 'I38.%' OR a.ICD10Code LIKE 'I39.%' THEN 'valvular disease'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I28.0', 'I28.8', 'I28.9') OR a.ICD10Code LIKE 'I26.%' OR a.ICD10Code LIKE 'I27.%' THEN 'pulmonary circulation disorders'
		
		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I73.1', 'I73.8', 'I73.9', 'I77.1', 'I79.0', 'I79.2', 'K55.1', 'K55.8', 'K55.9', 'Z95.8', 'Z95.9') OR a.ICD10Code LIKE 'I70.%' OR a.ICD10Code LIKE 'I71.%' THEN 'peripheral vascular disorders'
		
		WHEN a.ICD10Code LIKE 'I10.%' THEN 'hypertension, uncomplicated' 

		WHEN a.ICD10Code LIKE 'I11.%' OR a.ICD10Code LIKE 'I12.%' OR a.ICD10Code LIKE 'I13.%' OR a.ICD10Code LIKE 'I15.%' THEN 'hypertension, complicated'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('G04.1', 'G11.4', 'G80.1', 'G80.2', 'G83.0', 'G83.1', 'G83.2', 'G83.4', 'G83.9') OR a.ICD10Code LIKE 'G81.%' OR a.ICD10Code LIKE 'G82.%' THEN 'paralysis'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('G25.4', 'G25.5', 'G31.2', 'G31.8', 'G31.9', 'G93.1', 'G93.4', 'R47.0')
		OR SUBSTRING(a.ICD10Code, 1, 4) IN ('G10.', 'G11.', 'G12.', 'G13.', 'G20.', 'G21.', 'G22.', 'G32.', 'G35.' ,'G36.', 'G37.', 'G40.', 'G41.', 'R56.') THEN 'other neurological disorders' 

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I27.8', 'I27.9', 'J68.4', 'J70.1', 'J70.3')
		OR SUBSTRING(a.ICD10Code, 1, 4) IN ('J40.', 'J41.', 'J42.', 'J43.', 'J44.', 'J45.', 'J46.', 'J47.', 'J60.', 'J61.', 'J62.', 'J63.', 'J64.', 'J65.', 'J66.', 'J67.') THEN 'chronic pulmonary disease' 

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('E10.0', 'E10.1', 'E10.9', 'E11.0', 'E11.1', 'E11.9', 'E12.0', 'E12.1', 'E12.9', 'E13.0', 'E13.1', 'E13.9', 'E14.0', 'E14.1', 'E14.9') THEN 'diabetes, uncomplicated'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('E10.2', 'E10.8', 'E11.2', 'E11.3', 'E11.4', 'E11.5', 'E11.6', 'E11.7', 'E11.8', 'E12.2', 'E12.3', 'E12.4', 'E12.5', 'E12.6', 'E12.7', 'E12.8',
		'E13.2', 'E13.3', 'E13.4', 'E13.5', 'E13.6', 'E13.7', 'E13.8', 'E14.2', 'E14.3', 'E14.4', 'E14.5', 'E14.6', 'E14.7', 'E14.8') THEN 'diabetes, complicated'

		WHEN a.ICD10Code LIKE 'E00.%' OR a.ICD10Code LIKE 'E01.%' OR a.ICD10Code LIKE 'E02.%' OR a.ICD10Code LIKE 'E03.%' OR a.ICD10Code LIKE 'E89.0%' THEN 'hypothyroidism'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I12.0', 'I13.1', 'N25.0', 'Z49.0', 'Z49.1', 'Z49.2', 'Z94.0', 'Z99.2') OR a.ICD10Code LIKE 'N18.%' OR a.ICD10Code LIKE 'N19.%' THEN 'renal failure'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('I86.4', 'I98.2', 'K71.1', 'K71.3', 'K71.4', 'K71.5', 'K71.7', 'K76.0', 'K76.2', 'K76.3', 'K76.4', 'K76.5', 'K76.6', 'K76.7', 'K76.8', 'K76.9', 'Z94.4')
		OR a.ICD10Code LIKE 'B18.%' OR a.ICD10Code LIKE 'I85.%' OR a.ICD10Code LIKE 'K70.%' OR a.ICD10Code LIKE 'K72.%' OR a.ICD10Code LIKE 'K73.%' OR a.ICD10Code LIKE 'K74.%' THEN 'liver disease'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('K25.7', 'K25.9', 'K26.7', 'K26.9', 'K27.7', 'K27.9', 'K28.7', 'K28.9') THEN 'peptic ulcer disease'

		WHEN a.ICD10Code LIKE 'B20.%' OR a.ICD10Code LIKE 'B21.%' OR a.ICD10Code LIKE 'B22.%' OR a.ICD10Code LIKE 'B24.%' THEN 'aids/hiv'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('C90.0', 'C90.2') OR a.ICD10Code LIKE 'C81.%' OR a.ICD10Code LIKE 'C82.%' 
		OR a.ICD10Code LIKE 'C83.%' OR a.ICD10Code LIKE 'C84.%' OR a.ICD10Code LIKE 'C85.%' OR a.ICD10Code LIKE 'C88.%' OR a.ICD10Code LIKE 'C96.%' THEN 'lymphoma'

		WHEN SUBSTRING(a.ICD10Code, 1, 4) IN ('C77.', 'C78.', 'C79.', 'C80.') THEN 'metastatic cancer'

		WHEN SUBSTRING(a.ICD10Code, 1, 2) IN ('C0', 'C1', 'C6')
		OR SUBSTRING(a.ICD10Code, 1, 3) IN ('C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C30', 'C31', 'C32', 'C33', 'C34', 'C37', 'C48', 'C39', 
		'C40', 'C41', 'C42', 'C43', 'C45', 'C46', 'C47', 'C48', 'C49', 'C50', 'C51', 'C52', 'C53', 'C54', 'C55', 'C56', 'C57', 'C58', 'C70',
		'C71', 'C72', 'C73', 'C74', 'C75', 'C76', 'C97') THEN 'solid tumour without metastasis'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('L94.0', 'L94.1', 'L94.3', 'M12.0', 'M12.3', 'M31.0', 'M31.1', 'M31.2', 'M31.3', 'M46.1', 'M46.8', 'M46.9') OR
		SUBSTRING(a.ICD10Code, 1, 4) IN ('M05.', 'M06.', 'M08.', 'M30.', 'M32.', 'M33.', 'M34.', 'M35.', 'M45.') THEN 'rheumatoid arthritis'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('D69.1', 'D69.3', 'D69.4', 'D69.5', 'D69.6') OR SUBSTRING(a.ICD10Code, 1, 3) IN ('D65', 'D66', 'D67', 'D68') THEN 'coagulopathy'
 
		WHEN a.ICD10Code LIKE 'E66%' THEN 'obesity'

		WHEN SUBSTRING(a.ICD10Code, 1, 5)='R63.4' OR SUBSTRING(a.ICD10Code, 1, 3) IN ('E40', 'E41', 'E42', 'E43', 'E44', 'E45', 'E46', 'R64') THEN 'weight loss'

		WHEN SUBSTRING(a.ICD10Code, 1, 5)='E22.2' OR SUBSTRING(a.ICD10Code, 1, 3) IN ('E86', 'E87') THEN 'fluid and electrolyte disorders'

		WHEN SUBSTRING(a.ICD10Code, 1, 5)='D50.0' THEN 'blood loss anaemia'

		WHEN SUBSTRING(a.ICD10Code, 1, 5) IN ('D50.8', 'D50.9') OR SUBSTRING(a.ICD10Code, 1, 3) IN ('D51', 'D52', 'D53') THEN 'deficiency anaemia'

		ELSE NULL END AS type
	FROM CDWWork.Dim.ICD10 AS a
	LEFT OUTER JOIN OMHSP_PERC_Share.DOEx.Lookup_ICD10 AS lookup
		ON a.ICD10SID=lookup.ICD10SID
) AS z
WHERE type IS NOT NULL;


-- 3 minutes

SELECT
	PatientICN, RelativeTiming, type
INTO #relevant_dx
FROM(
	SELECT 
		#cohort.PatientICN, CASE WHEN diag.VisitDateTime<#cohort.VisitDateTime THEN 'lag1' ELSE 'cum3' END AS RelativeTiming, type 
	FROM CDWWork.Outpat.VDiagnosis AS diag
	INNER JOIN #ICD9
		ON diag.ICD9SID=#ICD9.ICD9SID
	INNER JOIN CDWWork.SPatient.SPatient AS pat
		ON diag.PatientSID=pat.PatientSID
	INNER JOIN #cohort
		ON pat.PatientICN=#cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, #cohort.VisitDateTime) AND DATEADD(day, 365*3, #cohort.VisitDateTime)
	WHERE diag.VisitDateTime>=CAST('2003-04-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2015-10-01 00:00:00' AS datetime2(0)) AND diag.ICD9SID>0
	UNION(
	SELECT 
		#cohort.PatientICN, CASE WHEN diag.VisitDateTime<#cohort.VisitDateTime THEN 'lag1' ELSE 'cum3' END AS RelativeTiming, type 
	FROM CDWWork.Outpat.VDiagnosis AS diag
	INNER JOIN #ICD10
		ON diag.ICD10SID=#ICD10.ICD10SID
	INNER JOIN CDWWork.SPatient.SPatient AS pat
		ON diag.PatientSID=pat.PatientSID
	INNER JOIN #cohort
		ON pat.PatientICN=#cohort.PatientICN AND diag.VisitDateTime BETWEEN DATEADD(day, -365, #cohort.VisitDateTime) AND DATEADD(day, 365*3, #cohort.VisitDateTime)
	WHERE diag.VisitDateTime>=CAST('2015-10-01 00:00:00' AS datetime2(0)) AND diag.VisitDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0)) AND diag.ICD10SID>0
	)
	UNION(
	SELECT
		#cohort.PatientICN, CASE WHEN diag.DischargeDateTime<#cohort.VisitDateTime THEN 'lag1' ELSE 'cum3' END AS RelativeTiming, type 
	FROM CDWWork.Inpat.InpatientDiagnosis AS diag
	INNER JOIN #ICD9
		ON diag.ICD9SID=#ICD9.ICD9SID
	INNER JOIN CDWWork.SPatient.SPatient AS pat
		ON diag.PatientSID=pat.PatientSID
	INNER JOIN #cohort
		ON pat.PatientICN=#cohort.PatientICN AND diag.DischargeDateTime BETWEEN DATEADD(day, -365, #cohort.VisitDateTime) AND DATEADD(day, 365*3, #cohort.VisitDateTime)
	WHERE diag.DischargeDateTime>=CAST('2003-04-01 00:00:00' AS datetime2(0)) AND diag.DischargeDateTime<CAST('2015-10-01 00:00:00' AS datetime2(0)) AND diag.ICD9SID>0
	)
	UNION(
	SELECT
		#cohort.PatientICN, CASE WHEN diag.DischargeDateTime<#cohort.VisitDateTime THEN 'lag1' ELSE 'cum3' END AS RelativeTiming, type 
	FROM CDWWork.Inpat.InpatientDiagnosis AS diag
	INNER JOIN #ICD10
		ON diag.ICD10SID=#ICD10.ICD10SID
	INNER JOIN CDWWork.SPatient.SPatient AS pat
		ON diag.PatientSID=pat.PatientSID
	INNER JOIN #cohort
		ON pat.PatientICN=#cohort.PatientICN AND diag.DischargeDateTime BETWEEN DATEADD(day, -365, #cohort.VisitDateTime) AND DATEADD(day, 365*3, #cohort.VisitDateTime)
	WHERE diag.DischargeDateTime>=CAST('2015-10-01 00:00:00' AS datetime2(0)) AND diag.DischargeDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0)) AND diag.ICD10SID>0
	)
) AS z


-- Elixhauser (1Y pre period only)
-- <1min
SELECT 
	PatientICN,
	9*congestive_heart_failure+6*pulmonary_circulation_disorders+3*peripheral_vascular_disorders
	-hypertension+5*paralysis+5*other_neurological_disorders+3*chronic_pulmonary_disease-3*diabetes_complicated
	+4*hypothyroidism+6*renal_failure+4*liver_disease+6*lymphoma+14*metastatic_cancer+7*solid_tumour
	+11*coagulopathy-5*obesity+9*weight_loss+11*fluid_and_electrolyte_disorders-3*blood_loss_anaemia-2*deficiency_anaemia
	-alcohol_abuse-5*psychoses-5*depression
	AS Elixhauser_AHRQ,
	7*congestive_heart_failure+5*cardiac_arrhythmia-valvular_disease+4*pulmonary_circulation_disorders+2*peripheral_vascular_disorders
	+7*paralysis+6*other_neurological_disorders+3*chronic_pulmonary_disease+5*renal_failure+11*liver_disease++9*lymphoma
	+12*metastatic_cancer+4*solid_tumour+3*coagulopathy-4*obesity+6*weight_loss+5*fluid_and_electrolyte_disorders
	-2*blood_loss_anaemia-2*deficiency_anaemia-3*depression
	AS Elixhauser_vanWalraven
FROM (
	SELECT
		PatientICN,
		MAX(CASE WHEN type='congestive heart failure' THEN 1 ELSE 0 END) AS congestive_heart_failure,
		MAX(CASE WHEN type='cardiac arrhythmia' THEN 1 ELSE 0 END) AS cardiac_arrhythmia,
		MAX(CASE WHEN type='valvular disease' THEN 1 ELSE 0 END) AS valvular_disease,
		MAX(CASE WHEN type='pulmonary circulation disorders' THEN 1 ELSE 0 END) AS pulmonary_circulation_disorders,
		MAX(CASE WHEN type='peripheral vascular disorders' THEN 1 ELSE 0 END) AS peripheral_vascular_disorders,
		MAX(CASE WHEN type IN ('hypertension, uncomplicated', 'hypertension, complicated') THEN 1 ELSE 0 END) AS hypertension,
		MAX(CASE WHEN type='paralysis' THEN 1 ELSE 0 END) AS paralysis,
		MAX(CASE WHEN type='other neurological disorders' THEN 1 ELSE 0 END) AS other_neurological_disorders,
		MAX(CASE WHEN type='chronic pulmonary disease' THEN 1 ELSE 0 END) AS chronic_pulmonary_disease,
		MAX(CASE WHEN type='diabetes, uncomplicated' THEN 1 ELSE 0 END) AS diabetes_uncomplicated,
		MAX(CASE WHEN type='diabetes, complicated' THEN 1 ELSE 0 END) AS diabetes_complicated,
		MAX(CASE WHEN type='hypothyroidism' THEN 1 ELSE 0 END) AS hypothyroidism,
		MAX(CASE WHEN type='renal failure' THEN 1 ELSE 0 END) AS renal_failure,
		MAX(CASE WHEN type='liver disease' THEN 1 ELSE 0 END) AS liver_disease,
		MAX(CASE WHEN type='peptic ulcer disease' THEN 1 ELSE 0 END) AS peptic_ulcer_disease,
		MAX(CASE WHEN type='aids/hiv' THEN 1 ELSE 0 END) AS aids_hiv,
		MAX(CASE WHEN type='lymphoma' THEN 1 ELSE 0 END) AS lymphoma,
		MAX(CASE WHEN type='metastatic cancer' THEN 1 ELSE 0 END) AS metastatic_cancer,
		MAX(CASE WHEN type='solid tumour without metastasis' THEN 1 ELSE 0 END) AS solid_tumour,
		MAX(CASE WHEN type='rheumatoid arthritis' THEN 1 ELSE 0 END) AS rheumatoid_arthritis,
		MAX(CASE WHEN type='coagulopathy' THEN 1 ELSE 0 END) AS coagulopathy,
		MAX(CASE WHEN type='obesity' THEN 1 ELSE 0 END) AS obesity,
		MAX(CASE WHEN type='weight loss' THEN 1 ELSE 0 END) AS weight_loss,
		MAX(CASE WHEN type='fluid and electrolyte disorders' THEN 1 ELSE 0 END) AS fluid_and_electrolyte_disorders,
		MAX(CASE WHEN type='blood loss anaemia' THEN 1 ELSE 0 END) AS blood_loss_anaemia,
		MAX(CASE WHEN type='deficiency anaemia' THEN 1 ELSE 0 END) AS deficiency_anaemia,
		MAX(CASE WHEN type='alcohol abuse' THEN 1 ELSE 0 END) AS alcohol_abuse,
		MAX(CASE WHEN type='psychoses' THEN 1 ELSE 0 END) AS psychoses,
		MAX(CASE WHEN type='depression' THEN 1 ELSE 0 END) AS depression
	FROM #relevant_dx
	WHERE RelativeTiming='lag1'
	GROUP BY PatientICN
) AS z")


healthoutcomes <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 
-- Diagnosis outcomes
-- includes CMS
-- <1min
SELECT
	PatientICN,
	MAX(CASE WHEN type='opioid overdose' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS opioid_overdose_lag1, MAX(CASE WHEN type='opioid overdose' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS opioid_overdose_cum3,
	MAX(CASE WHEN type='heroin overdose' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS heroin_overdose_lag1, MAX(CASE WHEN type='heroin overdose' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS heroin_overdose_cum3,
	MAX(CASE WHEN type='opioid dependence' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS OUD_lag1, MAX(CASE WHEN type='opioid dependence' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS OUD_cum3,
	MAX(CASE WHEN type='SUD' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS SUD_lag1, MAX(CASE WHEN type='SUD' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS SUD_cum3,
	MAX(CASE WHEN type='suicide' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS suicide_lag1, MAX(CASE WHEN type='suicide' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS suicide_cum3,
	MAX(CASE WHEN type='depression' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS depression_lag1, MAX(CASE WHEN type='depression' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS depression_cum3,
	MAX(CASE WHEN type='homeless' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS homeless_lag1, MAX(CASE WHEN type='homeless' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS homeless_cum3,
	MAX(CASE WHEN type='fall' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS fall_lag1, MAX(CASE WHEN type='fall' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS fall_cum3,
	MAX(CASE WHEN type='SAE' AND RelativeTiming='lag1' THEN 1 ELSE 0 END) AS SAE_lag1, MAX(CASE WHEN type='SAE' AND RelativeTiming='cum3' THEN 1 ELSE 0 END) AS SAE_cum3
FROM #relevant_dx
GROUP BY PatientICN")


pcp_claims <- merge(x=pcp_claims, y=elixhauser, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims <- merge(x=pcp_claims, y=healthoutcomes, by.x="PatientICN", by.y="PatientICN", all.x=T)

setnafill(pcp_claims, cols=c("Elixhauser_AHRQ", "Elixhauser_vanWalraven", 
                             "opioid_overdose_lag1", "opioid_overdose_cum3", "heroin_overdose_lag1", "heroin_overdose_cum3", "OUD_lag1", "OUD_cum3", "SUD_lag1", "SUD_cum3", "suicide_lag1", "suicide_cum3", "depression_lag1", "depression_cum3",
                             "homeless_lag1", "homeless_cum3", "fall_lag1", "fall_cum3", "SAE_lag1", "SAE_cum3"), fill=0)



elixhauser <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 

SELECT
PatientICN,
MAX(CASE WHEN type='congestive heart failure' THEN 1 ELSE 0 END) AS congestive_heart_failure,
MAX(CASE WHEN type='cardiac arrhythmia' THEN 1 ELSE 0 END) AS cardiac_arrhythmia,
MAX(CASE WHEN type='valvular disease' THEN 1 ELSE 0 END) AS valvular_disease,
MAX(CASE WHEN type='pulmonary circulation disorders' THEN 1 ELSE 0 END) AS pulmonary_circulation_disorders,
MAX(CASE WHEN type='peripheral vascular disorders' THEN 1 ELSE 0 END) AS peripheral_vascular_disorders,
MAX(CASE WHEN type IN ('hypertension, uncomplicated', 'hypertension, complicated') THEN 1 ELSE 0 END) AS hypertension,
MAX(CASE WHEN type='paralysis' THEN 1 ELSE 0 END) AS paralysis,
MAX(CASE WHEN type='other neurological disorders' THEN 1 ELSE 0 END) AS other_neurological_disorders,
MAX(CASE WHEN type='chronic pulmonary disease' THEN 1 ELSE 0 END) AS chronic_pulmonary_disease,
MAX(CASE WHEN type='diabetes, uncomplicated' THEN 1 ELSE 0 END) AS diabetes_uncomplicated,
MAX(CASE WHEN type='diabetes, complicated' THEN 1 ELSE 0 END) AS diabetes_complicated,
MAX(CASE WHEN type='hypothyroidism' THEN 1 ELSE 0 END) AS hypothyroidism,
MAX(CASE WHEN type='renal failure' THEN 1 ELSE 0 END) AS renal_failure,
MAX(CASE WHEN type='liver disease' THEN 1 ELSE 0 END) AS liver_disease,
MAX(CASE WHEN type='peptic ulcer disease' THEN 1 ELSE 0 END) AS peptic_ulcer_disease,
MAX(CASE WHEN type='aids/hiv' THEN 1 ELSE 0 END) AS aids_hiv,
MAX(CASE WHEN type='lymphoma' THEN 1 ELSE 0 END) AS lymphoma,
MAX(CASE WHEN type='metastatic cancer' THEN 1 ELSE 0 END) AS metastatic_cancer,
MAX(CASE WHEN type='solid tumour without metastasis' THEN 1 ELSE 0 END) AS solid_tumour,
MAX(CASE WHEN type='rheumatoid arthritis' THEN 1 ELSE 0 END) AS rheumatoid_arthritis,
MAX(CASE WHEN type='coagulopathy' THEN 1 ELSE 0 END) AS coagulopathy,
MAX(CASE WHEN type='obesity' THEN 1 ELSE 0 END) AS obesity,
MAX(CASE WHEN type='weight loss' THEN 1 ELSE 0 END) AS weight_loss,
MAX(CASE WHEN type='fluid and electrolyte disorders' THEN 1 ELSE 0 END) AS fluid_and_electrolyte_disorders,
MAX(CASE WHEN type='blood loss anaemia' THEN 1 ELSE 0 END) AS blood_loss_anaemia,
MAX(CASE WHEN type='deficiency anaemia' THEN 1 ELSE 0 END) AS deficiency_anaemia,
MAX(CASE WHEN type='alcohol abuse' THEN 1 ELSE 0 END) AS alcohol_abuse,
MAX(CASE WHEN type='psychoses' THEN 1 ELSE 0 END) AS psychoses,
MAX(CASE WHEN type='depression' THEN 1 ELSE 0 END) AS depression
FROM #relevant_dx
WHERE RelativeTiming='lag1'
GROUP BY PatientICN")

pcp_claims <- merge(x=pcp_claims, y=elixhauser, by.x="PatientICN", by.y="PatientICN", all.x=T)
setnafill(pcp_claims, cols=names(elixhauser), fill=0)

### BALANCE TABLE FIGURE 

pcp_claims$PeriodOfService <- as.factor(pcp_claims$PeriodOfService)
pcp_claims$PeriodOfService <- relevel(pcp_claims$PeriodOfService, "Gulf")
pcp_claims$MDC <- as.factor(pcp_claims$MDC)
pcp_claims$MDC <- relevel(pcp_claims$MDC, "Diseases of the circulatory system")
pcp_claims$RaceEthnicity <- as.factor(pcp_claims$RaceEthnicity)
pcp_claims$RaceEthnicity <- relevel(pcp_claims$RaceEthnicity, "White")
pcp_claims <- pcp_claims[!(EnrollmentPriority=="Other" | Married %in% c("", "4") | PeriodOfService=="" | MDC %in% c("Complications of pregnancy; childbirth; and the puerperium", "Congenital anomalies"))]

balance_MH <- felm(FE_MH_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+
                     congestive_heart_failure+cardiac_arrhythmia+valvular_disease+pulmonary_circulation_disorders+peripheral_vascular_disorders+hypertension+paralysis+other_neurological_disorders+chronic_pulmonary_disease+diabetes_uncomplicated+
                     diabetes_complicated+hypothyroidism+renal_failure+liver_disease+peptic_ulcer_disease+aids_hiv+lymphoma+metastatic_cancer+solid_tumour+rheumatoid_arthritis+coagulopathy+obesity+weight_loss+fluid_and_electrolyte_disorders+
                     blood_loss_anaemia+deficiency_anaemia+MH_dummy_lag1_v3+PH_dummy_lag1_v3+ACSC_lag1+WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

balance_PH <- felm(FE_PH_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+
                     congestive_heart_failure+cardiac_arrhythmia+valvular_disease+pulmonary_circulation_disorders+peripheral_vascular_disorders+hypertension+paralysis+other_neurological_disorders+chronic_pulmonary_disease+diabetes_uncomplicated+
                     diabetes_complicated+hypothyroidism+renal_failure+liver_disease+peptic_ulcer_disease+aids_hiv+lymphoma+metastatic_cancer+solid_tumour+rheumatoid_arthritis+coagulopathy+obesity+weight_loss+fluid_and_electrolyte_disorders+
                     blood_loss_anaemia+deficiency_anaemia+MH_dummy_lag1_v3+PH_dummy_lag1_v3+ACSC_lag1+WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

balance_ACSC <- felm(FE_ACSC_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+
                       congestive_heart_failure+cardiac_arrhythmia+valvular_disease+pulmonary_circulation_disorders+peripheral_vascular_disorders+hypertension+paralysis+other_neurological_disorders+chronic_pulmonary_disease+diabetes_uncomplicated+
                       diabetes_complicated+hypothyroidism+renal_failure+liver_disease+peptic_ulcer_disease+aids_hiv+lymphoma+metastatic_cancer+solid_tumour+rheumatoid_arthritis+coagulopathy+obesity+weight_loss+fluid_and_electrolyte_disorders+
                       blood_loss_anaemia+deficiency_anaemia+MH_dummy_lag1_v3+PH_dummy_lag1_v3+ACSC_lag1+WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

coef_names <- c("Age", "Race: Asian Pacific Islander", "Race: Black", "Race: Hispanic", "Race: Native American", "Race: Unknown", "Previously Married", "Never Married", "Log(1+Income)", "Medicare", "Medicaid", "Priority Group 2", "Priority Group 3", "Priority Group 4", "Priority Group 5", "Priority Group 6", "Priority Group 7", "Priority Group 8",
                "Period of Service: Korean War", "Period of Service: Other", "Period of Service: Post-Korean", "Period of Service: Post-Vietnam", "Period of Service: Vietnam", "Period of Service: WWII",
                "Service Connected Disability", "Unemployable", "Exposed to Agent Orange", "Exposed to Radiation", "Annual VA Check Amount", "Prior Year VHA Use", 
                "Prior Congestive Heart Failure", "Prior Cardiac Arrhythmia", "Prior Valcular Disease", "Prior Pulmonary Circulation Dis.", "Prior Peripheral Vascular Dis.", "Prior Hypertension", "Prior Paralysis", "Prior Other Neurological Dis.", "Prior Chronic Pulmonary", "Prior Diabetes (uncomp.)", "Prior Diabetes (comp.)",
                "Prior Hypothyroidism", "Prior Renal Failure", "Prior Liver Disease", "Prior Peptic Ulcer Disease", "Prior AIDS/HIV", "Prior Lymphoma", "Prior Metastatic Cancer", "Prior Solid Tumor", "Prior Rheumatoid Arthritis", "Prior Coagulopathy", 
                "Prior Obesity", "Prior Weight Loss", "Prior Fluid and Electrolyte Dis.", "Prior Blood Loss Anaemia", "Prior Deficiency Anaemia", 
                "Prior Year MH ED/Hosp/Psyc", "Prior Year PH ED/Hosp", "Prior ACSC Hosp", "Wait Time (days)", "MDC: Sym. & Signs & Ill-defined", "MDC: Blood & Blood-Forming", "MDC:Digestive", "MDC: Genitourinary",
                "MDC: Musc. & Connect. Tissue", "MDC: Nerv. Sys. & Sense Org.", "MDC: Respiratory", "MDC: Skin & Subcutaneous Tis.", "MDC: Endocr., Nutr. & Metab.",
                "MDC: Infect. & Paras.", "MDC: Injury & Pois.", "MDC: Mental Illness", "MDC: Neoplasms", "MDC: Residual codes, E codes")

layout(matrix(c(rep(1,(42)), rep(2,29), rep(3,29)), nrow=1, ncol=100, byrow=TRUE)) #save as 12 x 20

par(oma=c(1,0,0,0), mar=c(5.5,18,4.5,0.5))
plot(x=summary(balance_MH)$coef[,1], y=length(coef(balance_MH)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_MH))), pch=15, cex=1.25, col="royalblue4", yaxt="n", xaxt="n", cex.axis=1.8, xlab="", ylab="")
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_MH)$coef[,1]-1.96*summary(balance_MH)$coef[,2], length(coef(balance_MH)):1,
       summary(balance_MH)$coef[,1]+1.96*summary(balance_MH)$coef[,2], length(coef(balance_MH)):1, length=0.05, angle=90, code=3, col="royalblue4", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_MH)):1, x=-0.45, coef_names, adj=1, xpd=TRUE, cex=1.4)
text(y=length(coef(balance_MH)), x=0.25, paste("F-Stat: ", round(summary(balance_MH)$P.fstat[5],1)), cex=2, col="royalblue4")
title(main=list("Mental Health", cex=2.2, font=1, col="royalblue4"))

par(mar=c(5.5,0.5,4.5,0.5))
plot(x=summary(balance_PH)$coef[,1], y=length(coef(balance_PH)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_PH))), pch=15, cex=1.25, col="firebrick3", yaxt="n", xaxt="n", cex.axis=1.8, xlab="",)
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_PH)$coef[,1]-1.96*summary(balance_PH)$coef[,2], length(coef(balance_PH)):1,
       summary(balance_PH)$coef[,1]+1.96*summary(balance_PH)$coef[,2], length(coef(balance_PH)):1, length=0.05, angle=90, code=3, col="firebrick3", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_PH)), x=0.25, paste("F-Stat: ", round(summary(balance_PH)$P.fstat[5],1)), cex=2, col="firebrick3")
title(main=list("Circulatory Conditions", cex=2.2, font=1, col="firebrick3"))

par(mar=c(5.5,0.5,4.5,0.5))
plot(x=summary(balance_ACSC)$coef[,1], y=length(coef(balance_ACSC)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_ACSC))), pch=15, cex=1.25, col="darkorange", yaxt="n", xaxt="n", cex.axis=1.8, xlab="",)
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_ACSC)$coef[,1]-1.96*summary(balance_ACSC)$coef[,2], length(coef(balance_ACSC)):1,
       summary(balance_ACSC)$coef[,1]+1.96*summary(balance_ACSC)$coef[,2], length(coef(balance_ACSC)):1, length=0.05, angle=90, code=3, col="darkorange", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_ACSC)), x=0.25, paste("F-Stat: ", round(summary(balance_ACSC)$P.fstat[5],1)), cex=2.2, col="darkorange")
title(main=list("ACSC", cex=2, font=1, col="darkorange"))

mtext("Coefficient Estimates of Multivariate Regression", side=1, line=-1, outer=TRUE, cex=1.5, adj=0.61)

##### Chetty et al. (2014) figure:
# Step 1: Residualize all variables with baseline FEs
pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

resid_fn <- function(y){return(resid(feols(formula(paste(y, "~0|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins", sep="")), data=pcp_claims)))}

pcp_claims$Death3Year_resid <- resid_fn("Death3Year")
pcp_claims$herc3_resid <- resid_fn("herc3")
pcp_claims$MH_dummy_cum3_v1_resid <- resid_fn("MH_dummy_cum3_v1")
pcp_claims$PH_dummy_cum3_v1_resid <- resid_fn("PH_dummy_cum3_v1")
pcp_claims$ACSC_cum3_resid <- resid_fn("ACSC_cum3")

pcp_claims$Age_resid <- resid_fn("Age")
pcp_claims$Medicare_resid <- resid_fn("Medicare")
pcp_claims$Medicaid_resid <- resid_fn("Medicaid")
pcp_claims$ServiceConnectedFlag_resid <- resid_fn("ServiceConnectedFlag")
pcp_claims$UnemployableFlag_resid <- resid_fn("UnemployableFlag")
pcp_claims$ExposedToAgentOrangeFlag_resid <- resid_fn("ExposedToAgentOrangeFlag")
pcp_claims$RadiationExposureIndicatedFlag_resid <- resid_fn("RadiationExposureIndicatedFlag")
pcp_claims$AnnualVACheckAmount_resid <- resid_fn("AnnualVACheckAmount")

pcp_claims$EnrollmentPriority2 <- (pcp_claims$EnrollmentPriority=="GROUP 2"); pcp_claims$EnrollmentPriority3 <- (pcp_claims$EnrollmentPriority=="GROUP 3"); pcp_claims$EnrollmentPriority4 <- (pcp_claims$EnrollmentPriority=="GROUP 4"); 
pcp_claims$EnrollmentPriority5 <- (pcp_claims$EnrollmentPriority=="GROUP 5"); pcp_claims$EnrollmentPriority6 <- (pcp_claims$EnrollmentPriority=="GROUP 6"); pcp_claims$EnrollmentPriority7 <- (pcp_claims$EnrollmentPriority=="GROUP 7")
pcp_claims$EnrollmentPriority8 <- (pcp_claims$EnrollmentPriority=="GROUP 8");

pcp_claims$EnrollmentPriority2_resid <- resid_fn("EnrollmentPriority2")
pcp_claims$EnrollmentPriority3_resid <- resid_fn("EnrollmentPriority3")
pcp_claims$EnrollmentPriority4_resid <- resid_fn("EnrollmentPriority4")
pcp_claims$EnrollmentPriority5_resid <- resid_fn("EnrollmentPriority5")
pcp_claims$EnrollmentPriority6_resid <- resid_fn("EnrollmentPriority6")
pcp_claims$EnrollmentPriority7_resid <- resid_fn("EnrollmentPriority7")
pcp_claims$EnrollmentPriority8_resid <- resid_fn("EnrollmentPriority8")


pcp_claims$Elixhauser_resid <- resid_fn("Elixhauser_AHRQ")
pcp_claims$MH_dummy_lag1_v3_resid <- resid_fn("MH_dummy_lag1_v3")
pcp_claims$PH_dummy_lag1_v3_resid <- resid_fn("PH_dummy_lag1_v3")
pcp_claims$ACSC_lag1_resid <- resid_fn("ACSC_lag1")
pcp_claims$WaitTime_resid <- resid_fn("WaitTime")


# Step 2: Get predicted Y 
pcp_claims$MH_dummy_cum3_predicted <- felm(MH_dummy_cum3_v1_resid~Age_resid+Medicare_resid+Medicaid_resid+ServiceConnectedFlag_resid+UnemployableFlag_resid+ExposedToAgentOrangeFlag_resid+RadiationExposureIndicatedFlag_resid+AnnualVACheckAmount_resid
                                           +EnrollmentPriority2_resid+EnrollmentPriority3_resid+EnrollmentPriority4_resid+EnrollmentPriority5_resid+EnrollmentPriority6_resid+EnrollmentPriority7_resid+EnrollmentPriority8_resid
                                           +Elixhauser_resid+ACSC_lag1_resid+MH_dummy_lag1_v3_resid+PH_dummy_lag1_v3_resid+WaitTime_resid, data=pcp_claims)$fitted.values

pcp_claims$PH_dummy_cum3_predicted <- felm(PH_dummy_cum3_v1_resid~Age_resid+Medicare_resid+Medicaid_resid+ServiceConnectedFlag_resid+UnemployableFlag_resid+ExposedToAgentOrangeFlag_resid+RadiationExposureIndicatedFlag_resid+AnnualVACheckAmount_resid
                                           +EnrollmentPriority2_resid+EnrollmentPriority3_resid+EnrollmentPriority4_resid+EnrollmentPriority5_resid+EnrollmentPriority6_resid+EnrollmentPriority7_resid+EnrollmentPriority8_resid
                                           +Elixhauser_resid+ACSC_lag1_resid+MH_dummy_lag1_v3_resid+PH_dummy_lag1_v3_resid+WaitTime_resid, data=pcp_claims)$fitted.values

pcp_claims$ACSC_cum3_predicted <- felm(ACSC_cum3_resid~Age_resid+Medicare_resid+Medicaid_resid+ServiceConnectedFlag_resid+UnemployableFlag_resid+ExposedToAgentOrangeFlag_resid+RadiationExposureIndicatedFlag_resid+AnnualVACheckAmount_resid
                                       +EnrollmentPriority2_resid+EnrollmentPriority3_resid+EnrollmentPriority4_resid+EnrollmentPriority5_resid+EnrollmentPriority6_resid+EnrollmentPriority7_resid+EnrollmentPriority8_resid
                                       +Elixhauser_resid+ACSC_lag1_resid+MH_dummy_lag1_v3_resid+PH_dummy_lag1_v3_resid+WaitTime_resid, data=pcp_claims)$fitted.values

pcp_claims$herc3_predicted <- felm(herc3~Age_resid+Medicare_resid+Medicaid_resid+ServiceConnectedFlag_resid+UnemployableFlag_resid+ExposedToAgentOrangeFlag_resid+RadiationExposureIndicatedFlag_resid+AnnualVACheckAmount_resid
                                       +EnrollmentPriority2_resid+EnrollmentPriority3_resid+EnrollmentPriority4_resid+EnrollmentPriority5_resid+EnrollmentPriority6_resid+EnrollmentPriority7_resid+EnrollmentPriority8_resid
                                       +Elixhauser_resid+ACSC_lag1_resid+MH_dummy_lag1_v3_resid+PH_dummy_lag1_v3_resid+WaitTime_resid, data=pcp_claims)$fitted.values

pcp_claims$Death3Year_predicted <- felm(Death3Year_resid~Age_resid+Medicare_resid+Medicaid_resid+ServiceConnectedFlag_resid+UnemployableFlag_resid+ExposedToAgentOrangeFlag_resid+RadiationExposureIndicatedFlag_resid+AnnualVACheckAmount_resid
                                       +EnrollmentPriority2_resid+EnrollmentPriority3_resid+EnrollmentPriority4_resid+EnrollmentPriority5_resid+EnrollmentPriority6_resid+EnrollmentPriority7_resid+EnrollmentPriority8_resid
                                       +Elixhauser_resid+ACSC_lag1_resid+MH_dummy_lag1_v3_resid+PH_dummy_lag1_v3_resid+WaitTime_resid, data=pcp_claims)$fitted.values


# Step 3: Residualize Y and X (effectiveness)
pcp_claims$MH_dummy_cum3_predicted <- resid_fn("MH_dummy_cum3_predicted")
pcp_claims$PH_dummy_cum3_predicted <- resid_fn("PH_dummy_cum3_predicted")
pcp_claims$ACSC_cum3_predicted <- resid_fn("ACSC_cum3_predicted")
pcp_claims$Death3Year_predicted <- resid_fn("Death3Year_predicted")
pcp_claims$herc3_predicted <- resid_fn("herc3_predicted")
pcp_claims$FE_MH_resid <- -resid_fn("FE_MH")
pcp_claims$FE_PH_resid <- -resid_fn("FE_PH")
pcp_claims$FE_ACSC_resid <- -resid_fn("FE_ACSC")

# Step 4: Split X into 20 evenly sized bins
chetty <- pcp_claims[, c("FE_MH_resid", "FE_PH_resid", "FE_ACSC_resid", "MH_dummy_cum3_predicted", "MH_dummy_cum3_v1_resid", "PH_dummy_cum3_predicted", "PH_dummy_cum3_v1_resid", "ACSC_cum3_predicted", "ACSC_cum3_resid", 
                         "Death3Year_resid", "Death3Year_predicted", "herc3_predicted", "herc3_resid")]
chetty$FE_MH_bins <- cut2(chetty$FE_MH_resid, g=20); chetty$FE_PH_bins <- cut2(chetty$FE_PH_resid, g=20); chetty$FE_ACSC_bins <- cut2(chetty$FE_ACSC_resid, g=20)
chetty_MH <- chetty[, .(FE_MH_resid=mean(FE_MH_resid), MH_dummy_cum3_predicted=mean(MH_dummy_cum3_predicted), MH_dummy_cum3_v1_resid=mean(MH_dummy_cum3_v1_resid), 
                        Death3Year_predicted=mean(Death3Year_predicted), Death3Year_resid=mean(Death3Year_resid), herc3_predicted=mean(herc3_predicted), herc3_resid=mean(herc3_resid)), by=c("FE_MH_bins")]
chetty_PH <- chetty[, .(FE_PH_resid=mean(FE_PH_resid), PH_dummy_cum3_predicted=mean(PH_dummy_cum3_predicted), PH_dummy_cum3_v1_resid=mean(PH_dummy_cum3_v1_resid), 
                        Death3Year_predicted=mean(Death3Year_predicted), Death3Year_resid=mean(Death3Year_resid), herc3_predicted=mean(herc3_predicted), herc3_resid=mean(herc3_resid)), by=c("FE_PH_bins")]
chetty_ACSC <- chetty[, .(FE_ACSC_resid=mean(FE_ACSC_resid), ACSC_cum3_predicted=mean(ACSC_cum3_predicted), ACSC_cum3_resid=mean(ACSC_cum3_resid), 
                          Death3Year_predicted=mean(Death3Year_predicted), Death3Year_resid=mean(Death3Year_resid), herc3_predicted=mean(herc3_predicted), herc3_resid=mean(herc3_resid)), by=c("FE_ACSC_bins")]


# Step 5: Plot
par(mfrow=c(2,2), mar=c(4.5,4.7,3,2))
plot(chetty_MH[,c("FE_MH_resid", "MH_dummy_cum3_v1_resid")], cex=1.5, xlab="Mental Health Effectiveness", ylab="3 Year Mental Health ED/Hospitalization", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_MH[,c("FE_MH_resid", "MH_dummy_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=-0.03, y=-0.005, paste("Coef: ", round(summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=-0.03, y=0.01, col="red", paste("Coef: ", round(summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_PH[,c("FE_PH_resid", "PH_dummy_cum3_v1_resid")], cex=1.5, xlab="Circulatory Condition Effectiveness", ylab="3 Year Circulatory ED/Hospitalization", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_PH[,c("FE_PH_resid", "PH_dummy_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=-0.045, y=-0.005, paste("Coef: ", round(summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=-0.045, y=0.015, col="red", paste("Coef: ", round(summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_ACSC[,c("FE_ACSC_resid", "ACSC_cum3_resid")], cex=1.5, xlab="ACSC Effectiveness", ylab="3 Year ACSC Hospitalizations", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_ACSC[,c("FE_ACSC_resid", "ACSC_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=-0.032, y=-0.005, paste("Coef: ", round(summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=-0.032, y=0.01, col="red", paste("Coef: ", round(summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)



######################################################
############      Appendix Figure A2      ############
######################################################

balance_MH <- felm(FE_MH_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+
                     WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[PriorVisit==0])

balance_PH <- felm(FE_PH_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+
                     WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[PriorVisit==0])

balance_ACSC <- felm(FE_ACSC_bare_std~Age+RaceEthnicity+factor(Married)+log(1+PatientIncome)+Medicare+Medicaid+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+
                       WaitTime+MDC|InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[PriorVisit==0])

coef_names <- c("Age", "Race: Asian Pacific Islander", "Race: Black", "Race: Hispanic", "Race: Native American", "Race: Unknown", "Previously Married", "Never Married", "Log(1+Income)", "Medicare", "Medicaid", "Priority Group 2", "Priority Group 3", "Priority Group 4", "Priority Group 5", "Priority Group 6", "Priority Group 7", "Priority Group 8",
                "Period of Service: Korean War", "Period of Service: Other", "Period of Service: Post-Korean", "Period of Service: Post-Vietnam", "Period of Service: Vietnam", "Period of Service: WWII",
                "Service Connected Disability", "Unemployable", "Exposed to Agent Orange", "Exposed to Radiation", "Annual VA Check Amount", "Wait Time (days)", "MDC: Sym. & Signs & Ill-defined", "MDC: Blood & Blood-Forming", "MDC:Digestive", "MDC: Genitourinary",
                "MDC: Musc. & Connect. Tissue", "MDC: Nerv. Sys. & Sense Org.", "MDC: Respiratory", "MDC: Skin & Subcutaneous Tis.", "MDC: Endocr., Nutr. & Metab.",
                "MDC: Infect. & Paras.", "MDC: Injury & Pois.", "MDC: Mental Illness", "MDC: Neoplasms", "MDC: Residual codes, E codes")

layout(matrix(c(rep(1,(42)), rep(2,29), rep(3,29)), nrow=1, ncol=100, byrow=TRUE)) #save as 12 x 20

par(oma=c(1,0,0,0), mar=c(5.5,18,4.5,0.5))
plot(x=summary(balance_MH)$coef[,1], y=length(coef(balance_MH)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_MH))), pch=15, cex=1.25, col="royalblue4", yaxt="n", xaxt="n", cex.axis=1.8, xlab="", ylab="")
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_MH)$coef[,1]-1.96*summary(balance_MH)$coef[,2], length(coef(balance_MH)):1,
       summary(balance_MH)$coef[,1]+1.96*summary(balance_MH)$coef[,2], length(coef(balance_MH)):1, length=0.05, angle=90, code=3, col="royalblue4", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_MH)):1, x=-0.45, coef_names, adj=1, xpd=TRUE, cex=1.4)
text(y=length(coef(balance_MH)), x=0.25, paste("F-Stat: ", round(summary(balance_MH)$P.fstat[5],1)), cex=2, col="royalblue4")
title(main=list("Mental Health", cex=2.2, font=1, col="royalblue4"))

par(mar=c(5.5,0.5,4.5,0.5))
plot(x=summary(balance_PH)$coef[,1], y=length(coef(balance_PH)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_PH))), pch=15, cex=1.25, col="firebrick3", yaxt="n", xaxt="n", cex.axis=1.8, xlab="",)
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_PH)$coef[,1]-1.96*summary(balance_PH)$coef[,2], length(coef(balance_PH)):1,
       summary(balance_PH)$coef[,1]+1.96*summary(balance_PH)$coef[,2], length(coef(balance_PH)):1, length=0.05, angle=90, code=3, col="firebrick3", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_PH)), x=0.25, paste("F-Stat: ", round(summary(balance_PH)$P.fstat[5],1)), cex=2, col="firebrick3")
title(main=list("Circulatory Conditions", cex=2.2, font=1, col="firebrick3"))

par(mar=c(5.5,0.5,4.5,0.5))
plot(x=summary(balance_ACSC)$coef[,1], y=length(coef(balance_ACSC)):1, xlim=c(-0.4,0.4), ylim=c(0,length(coef(balance_ACSC))), pch=15, cex=1.25, col="darkorange", yaxt="n", xaxt="n", cex.axis=1.8, xlab="",)
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), c(-0.4, -0.2,0,0.2,0.4), cex.axis=2)
arrows(summary(balance_ACSC)$coef[,1]-1.96*summary(balance_ACSC)$coef[,2], length(coef(balance_ACSC)):1,
       summary(balance_ACSC)$coef[,1]+1.96*summary(balance_ACSC)$coef[,2], length(coef(balance_ACSC)):1, length=0.05, angle=90, code=3, col="darkorange", lwd=1.5)
abline(v=0, lwd=1.5, lty=2)
text(y=length(coef(balance_ACSC)), x=0.25, paste("F-Stat: ", round(summary(balance_ACSC)$P.fstat[5],1)), cex=2.2, col="darkorange")
title(main=list("ACSC", cex=2, font=1, col="darkorange"))

mtext("Coefficient Estimates of Multivariate Regression", side=1, line=-1, outer=TRUE, cex=1.5, adj=0.61)


# Chetty et al. (2014)
chetty <- pcp_claims[PriorVisit==0, c("FE_MH_resid", "FE_PH_resid", "FE_ACSC_resid", "MH_dummy_cum3_predicted", "MH_dummy_cum3_v1_resid", "PH_dummy_cum3_predicted", "PH_dummy_cum3_v1_resid", "ACSC_cum3_predicted", "ACSC_cum3_resid")]
chetty$FE_MH_bins <- cut2(chetty$FE_MH_resid, g=20); chetty$FE_PH_bins <- cut2(chetty$FE_PH_resid, g=20); chetty$FE_ACSC_bins <- cut2(chetty$FE_ACSC_resid, g=20)
chetty_MH <- chetty[, .(FE_MH_resid=mean(FE_MH_resid), MH_dummy_cum3_predicted=mean(MH_dummy_cum3_predicted), MH_dummy_cum3_v1_resid=mean(MH_dummy_cum3_v1_resid)), by=c("FE_MH_bins")]
chetty_PH <- chetty[, .(FE_PH_resid=mean(FE_PH_resid), PH_dummy_cum3_predicted=mean(PH_dummy_cum3_predicted), PH_dummy_cum3_v1_resid=mean(PH_dummy_cum3_v1_resid)), by=c("FE_PH_bins")]
chetty_ACSC <- chetty[, .(FE_ACSC_resid=mean(FE_ACSC_resid), ACSC_cum3_predicted=mean(ACSC_cum3_predicted), ACSC_cum3_resid=mean(ACSC_cum3_resid)), by=c("FE_ACSC_bins")]

# Plot
par(mfrow=c(2,2), mar=c(4.5,4.7,3,2))
plot(chetty_MH[,c("FE_MH_resid", "MH_dummy_cum3_v1_resid")], cex=1.5, xlab="Mental Health Effectiveness", ylab="3 Year Mental Health ED/Hospitalization", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_MH[,c("FE_MH_resid", "MH_dummy_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0][PriorVisit==0]))$coef[2], lwd=2)
text(x=-0.03, y=-0.005, paste("Coef: ", round(summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(MH_dummy_cum3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2], lwd=2, col="red")
text(x=-0.03, y=0.01, col="red", paste("Coef: ", round(summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(MH_dummy_cum3_v1_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_PH[,c("FE_PH_resid", "PH_dummy_cum3_v1_resid")], cex=1.5, xlab="Circulatory Condition Effectiveness", ylab="3 Year Circulatory ED/Hospitalization", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_PH[,c("FE_PH_resid", "PH_dummy_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2], lwd=2)
text(x=-0.045, y=-0.005, paste("Coef: ", round(summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(PH_dummy_cum3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2], lwd=2, col="red")
text(x=-0.045, y=0.01, col="red", paste("Coef: ", round(summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(PH_dummy_cum3_v1_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_ACSC[,c("FE_ACSC_resid", "ACSC_cum3_resid")], cex=1.5, xlab="ACSC Effectiveness", ylab="3 Year ACSC Hospitalizations", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_ACSC[,c("FE_ACSC_resid", "ACSC_cum3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2], lwd=2)
text(x=-0.032, y=-0.005, paste("Coef: ", round(summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(ACSC_cum3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[1], b=summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2], lwd=2, col="red")
text(x=-0.032, y=0.01, col="red", paste("Coef: ", round(summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2],3), " (", round(summary(felm(ACSC_cum3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims[PriorVisit==0]))$coef[2,2], 3), ")", sep=""), cex=1.3)




###########################################################
#############        APPENDIX FIGURE A3        ############
###########################################################

par(mfrow=c(2,2), mar=c(4.5,4.7,3,2))
plot(chetty_MH[,c("FE_MH_resid", "Death3Year_resid")], cex=1.5, xlab="Mental Health Effectiveness", ylab="3 Year Mortality", cex.lab=1.5, cex.axis=1.5, col="red", pch=17, ylim=c(-0.008,0.008))
points(chetty_MH[,c("FE_MH_resid", "Death3Year_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(Death3Year_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.02, y=0.001, paste("Coef: ", round(summary(felm(Death3Year_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(Death3Year_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.02, y=-0.005, col="red", paste("Coef: ", round(summary(felm(Death3Year_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_PH[,c("FE_PH_resid", "Death3Year_resid")], cex=1.5, xlab="Circulatory Condition Effectiveness", ylab="3 Year Mortality", cex.lab=1.5, cex.axis=1.5, col="red", pch=17, ylim=c(-0.01,0.01))
points(chetty_PH[,c("FE_PH_resid", "Death3Year_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(Death3Year_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.035, y=0.001, paste("Coef: ", round(summary(felm(Death3Year_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(Death3Year_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.035, y=-0.005, col="red", paste("Coef: ", round(summary(felm(Death3Year_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_ACSC[,c("FE_ACSC_resid", "Death3Year_resid")], cex=1.5, xlab="ACSC Effectiveness", ylab="3 Year Mortality", cex.lab=1.5, cex.axis=1.5, col="red", pch=17, ylim=c(-0.01,0.01))
points(chetty_ACSC[,c("FE_ACSC_resid", "Death3Year_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(Death3Year_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.02, y=0.001, paste("Coef: ", round(summary(felm(Death3Year_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(Death3Year_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(Death3Year_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.02, y=-0.005, col="red", paste("Coef: ", round(summary(felm(Death3Year_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(Death3Year_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)



par(mfrow=c(2,2), mar=c(4.5,4.7,3,2))
plot(chetty_MH[,c("FE_MH_resid", "herc3_resid")], cex=1.5, xlab="Mental Health Effectiveness", ylab="3 Year Log Spending", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_MH[,c("FE_MH_resid", "herc3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(herc3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.02, y=0.01, paste("Coef: ", round(summary(felm(herc3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_predicted~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(herc3_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.02, y=-0.05, col="red", paste("Coef: ", round(summary(felm(herc3_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_resid~FE_MH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_PH[,c("FE_PH_resid", "herc3_resid")], cex=1.5, xlab="Circulatory Condition Effectiveness", ylab="3 Year Log Spending", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_PH[,c("FE_PH_resid", "herc3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(herc3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.035, y=0.01, paste("Coef: ", round(summary(felm(herc3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_predicted~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(herc3_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.035, y=-0.05, col="red", paste("Coef: ", round(summary(felm(herc3_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_resid~FE_PH_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)

plot(chetty_ACSC[,c("FE_ACSC_resid", "herc3_resid")], cex=1.5, xlab="ACSC Effectiveness", ylab="3 Year Log Spending", cex.lab=1.5, cex.axis=1.5, col="red", pch=17)
points(chetty_ACSC[,c("FE_ACSC_resid", "herc3_predicted")], pch=19, cex=1.5)
abline(a=summary(felm(herc3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2)
text(x=0.02, y=0.01, paste("Coef: ", round(summary(felm(herc3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_predicted~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)
abline(a=summary(felm(herc3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[1], b=summary(felm(herc3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2], lwd=2, col="red")
text(x=0.02, y=-0.05, col="red", paste("Coef: ", round(summary(felm(herc3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2],3), " (", round(summary(felm(herc3_resid~FE_ACSC_resid|0|0|InstitutionSID, data=pcp_claims))$coef[2,2], 3), ")", sep=""), cex=1.3)


